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Abstract 

The Method of Invariant Grid (MIG) is an iterative procedure for model reduction 
in chemical kinetics which is based on the notion of Slow Invariant Manifold (SIM) 
[l]-[4]. Important role, in that method, is played by the initial grid which, once re- 
fined, gives a description of the invariant manifold: the invariant grid. A convenient 
way to get a first approximation of the SIM is given by the Spectral Quasi Equi- 
librium Manifold (SQEM) [l]-[2]. In the present paper, a flexible numerical method 
to construct the discrete analog of a Quasi Equilibrium Manifold, in any dimen- 
sion, is presented. That object is named Quasi Equilibrium Grid (QEG), while the 
procedure Quasi Equilibrium Grid Algorithm. Extensions of the QEM notion are 
also suggested. The QEG is a numerical tool which can be used to find a grid-based 
approximation for the locus of minima of a convex function under some linear con- 
straints. The method is validated by construction of one and two-dimensional grids 
for model hydrogen oxidation reaction. 

Key words: Chemical kinetics, model reduction, invariant manifold, entropy, 
nonlinear dynamics, Lagrange multipliers method, variational problem. 



1 Introduction 



Relaxation of complex systems is often characterized by a fast dynamics dur- 
ing a short initial stage, while the remaining period lasts much longer and 
it evolves along low- dimensional surfaces in the phase space known as Slow 
Invariant Manifolds (SIM). In that scenario, a simplified macroscopic descrip- 
tion of a complex system can be attained by extracting only the slow dy- 
namics and neglecting the fast one. For this reason, much effort was spent 
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to develop model reduction methods (the Method of Invariant Grid (MIG) 
[T1I21I3B] . the Intrinsic Low Dimensional Manifold method (ILDM) [9lfT0] . the 
Computational Singular Perturbation method (CSP) [TTf 121113] . etc.) based on 
the notion of SIM. The introduction of a convex Lyapunov function G, when- 
ever the complex system is supported by such a function, also proves to be 
very helpful in model reduction [Hl8] . Indeed, it was shown that, through a G 
function, good approximations of the SIM can be found (e.g. by constructing 
the Spectral Quasi Equilibrium Manifold - SQEM - or the Symmetric En- 
tropic Intrinsic Low Dimensional Manifold - SEILDM - jTH2"] ) and refined by 
some efficient MIG iterations. Moreover, it has been shown that the notion 
of QEM is also very useful in different fields. For example, it was used in the 
implementation of Lattice Boltzmann schemes [517] . Construction of a QEM 
is analytically possible by using the Lagrange multipliers method. However, 
its implementation becomes too complicated as soon as the number of vari- 
ables of the problem becomes large: efficient methods, for constructing large 
dimensional QEM, are still missing. Therefore, in the present paper the no- 
tion of Quasi Equilibrium Grid (QEG) will be introduced, as a discrete analog 
of QEM, and a constructive algorithm, applicable in any dimension, will be 
developed. The procedure suggested proves to be a very flexible tool, so it 
is possible to get some other SIM approximations, all based on the previous 
algorithm, even more accurate than the QEG itself. 



2 Paper organization 

The paper is organized as follows. In Section [3], some basic notions are out- 
lined: in particular, the general equations of dissipative reaction kinetics are 
reviewed, in the notations which are used throughout the paper. At the end of 
that Section, the Method of Invariant Grid (MIG) and the notion of thermo- 
dynamic projector are briefly discussed (Section 13.21) . In Section HI the QEM 
definition and its geometrical interpretation is given, while in Section [5] the 
ID Quasi Equilibrium Grid Algorithm is presented. That algorithm is also il- 
lustrated, by means of an example, in Section The ID Algorithm extension 
to multi-dimensional grids is developed in Section [7J In particular, two pos- 
sible extension strategies are analyzed: the straightforward extension (Section 
I7.ip and, by following the general idea given in [3], the flag extension (Sec- 
tion I7.2p . Here, it is also shown how the flexibility of the flag extension allows 
to get SIM approximation which is better than the Spectral-QEG (Section 
IT. 3p : the notions of Guided- QEG and Symmetric Entropic Guided- QEG are 
introduced. An illustrative example, in Section [BJ shows how those different 
extension techniques work in practice. In order to find out how accurate is 
their SIM description, they are also compared on the base of the invariance 
defect (Section l8.3p . Finally, results are discussed in Section O 



2 



3 Theoretical background 

3.1 Dissipative reaction kinetics 

In a closed system with n chemical species A\, A n , participating in a com- 
plex reaction, a generic reversible reaction step can be written as a stoichio- 
metric equation: 

a s iAx + ... + a sn A n ^± (3 sX Ai + ... + f3 sn A n , (1) 

where s is the reaction index, s = 1, r (r steps in total), and the integers a si 
and (3 S i are stoichiometric coefficients of the step s. For each reaction step, we 
can introduce n-component vectors a s and (3 S , with components a S i and /3 S j, 
and the stoichiometric vector 7 s =/3 s -a s . For every A4 the extensive variable 
Ni describes the number of particles of that species. If V is the volume, then 
the concentration of A» is c, = Ni/V. Dynamics of the species concentration 
according to the stoichiometric mechanism (pQ) reads: 

N=VJ(c), J(c) = Z:=il s W s (c), (2) 

where dot denotes the time derivative and W s (c) is the reaction rate function 
of the step s. In particular, the polynomial form of the reaction rate function 
is provided by the mass action law: 

W s (c) = W+(c) - W-(c) = kt(T) f[c?-K(T) f[<* , (3) 

i=l i=l 

where kf(T) and k~(T) are the constants of the direct and of the inverse 
reactions rates of the step s respectively. The most popular form of their 
dependence is given by the Arrhenius equation: 

kf(T) = afT b * exp(S±/k B )exp(-H±/k B T). 

In the latter equation, af,bf are constants and Hf, Sf activation enthalpies 
and entropies respectively. The rate constants are not independent. Indeed, 
the principle of detail balance gives a relation between these quantities: 

W+(c e «) = W s -{c e o), Vs = 1, ...,r, (4) 

where the positive vector c eq (T) is the equilibrium of the system (j2J). In order 
to obtain a closed system of equations, one should supply an equation for 
the volume V. For an isolated system the extra-equations are U, V = const 
(where U is the internal energy), for an isochoric isothermal system we get 
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V, T = const, and so forth. For example, equation ([2]) in the latter case simply 
takes the form: 

c = ir ls W s {c) = J(c). (5) 

s=l 

Finally, also other linear constraints, related to the conservation of atoms, 
must be considered. In general such conservation laws can have the following 
form: 

Dc = const, (6) 

where I fixed and linearly independent vectors are the rows of the / x n 
matrix D, and const is a constant vector. 



3.2 Outline of the method of invariant grid 



In this section, we give an outline of the MIG for chemical kinetics. For details 
see Refs. 



3.2.1 Thermodynamic potential 

If we turn our attention to perfectly stirred closed chemically active mixtures, 
then dissipative properties of such systems can be characterized with a ther- 
modynamic potential which is the Lyapunov function of equation (121) . That 
function implements Second Law of thermodynamics: it means that during the 
concentrations evolution in time, from the initial condition to the equilibrium 
state, the Lyapunov function must decrease monotonically. Therefore if G(c) 
is the Lyapunov function, c eg (equilibrium state) is its point of global mini- 
mum in the phase space. A simple example of a function G is given by the free 
energy of ideal gas in a constant volume and under a constant temperature: 

G = f>[m( Ci /cf)-l]. (7) 

8=1 

When G is known, also its gradient VG and the matrix of second derivatives 
H =|| d 2 G I 'dcidcj || can be evaluated, so that it is possible to introduce the 
thermodynamic scalar product as follows: 

(x,y) = (x,Hy), (8) 

where the notation (, ) is the usual Euclidean scalar product. 
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3.2.2 The invariance condition 



Let us consider ft as a manifold of a reduced description. The invariance 
requirement reads: 

c(o) eO^ c(t) efi, vt > o. (9) 

Let P be a projector on the tangent bundle of the manifold ft. The manifold 
ft is invariant with respect to the system (jSJ) if and only if the following 
invariance equation (IE) holds: 

[1 - P] J(c) = 0, Vc G ft. (10) 

When the manifold is not invariant, it is not able to satisfy the invariance 
condition so that: 

3c : A = [1 - P]J(co) ± 0, (11) 

where Ao is the defect of invariance. One way to find the SIM is to solve the 
IE iteratively starting from an appropriate initial manifold. 



3.2.3 Thermodynamic 'projector 

Let us now discuss further the projector appearing in the invariance equation. 
It is an operator which for each point c G ft projects the vectors J(c) onto the 
tangent subspace of the manifold producing, in this way, the induced vector 
field PJ(c). In general, condition ffTUl) does not require any special constraint 
for the projector P. However, the thermodynamic properties of the kinetic 
equations (j2J) define the projector unambiguously [THUS] . To this end, let us 
define a differential of G, that is linear functional: 

DG(x) = (VG(c), x). (12) 

A special class of projectors is the thermodynamic one. If a projector belongs 
to this class then the induced vector field respects the dissipation inequality: 

DG(PJ) < 0, Vc G ft. (13) 
It has been shown that a projector P respects the f|T3|) if and only if [8J: 

ker P C ker DG, Vc G ft, (14) 

where ker denotes the null-space of an operator. It is clear now that if one 
wants to solve equation ([10]) . then a projector must be specified. Here we 
remind the way to construct the thermodynamic projector which will be used 
in MIG procedure pp. This projector depends on the concentration point c 
and on the tangent space to the manifold ft. 
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We are looking for a grid approximation of a q- dimensional SIM. Let Q be 
a discrete subset of a g- dimensional parameter space R q and let F \g be a 
mapping of Q into the concentration space. If we select an approximation 
procedure to restore the smooth map F from the discrete map F \g (we need 
a very small part of F, derivatives of F in the grid points only), then the 
derivatives f i = dF/dr/i are available, and for each grid point the tangent 
space is: 

T y = Lin{f t }, z = l,...,n. (15) 

We assume that one of points y G Q maps into the equilibrium, and in 
other points intersection of the manifold with G levels is transversal (i.e. 
(DG)F{ y ){x) 7^ for some x G T y ). Let us consider the subspace T 0y = 
(T y Piker DG). In order to define the thermodynamic projector, it is required, 
if T 0y ^ T y , to introduce the vector e y which satisfies the following conditions: 

e G T 

Cy C ±y, 

< (e y ,x) = 0, \/x G T 0y , 
DG(e y ) = 1. 

Let Po be the orthogonal projector on T 0y with respect to the entropic scalar 
product dSJ), then the thermodynamic projection of a vector x is defined as: 

T y ^T v ^Px = P x + e y DG{x) 

(16) 

T y =Ty^PX = P Q X. 



3.2.4 Iterative procedures: the Newton method with incomplete linearization 

When MIG method is applied, not a manifold is searched as a solution, but a 
set of concentration points whose defect of invariance is sufficiently small: let f2 
denote that solution (invariant grid). MIG is an iterative procedure: this means 
that, at the beginning, only an initial approximation f2 of f2 is available. In 
general, Qq does not respect the invariance condition ffTUl) satisfactorily so the 
pip holds: for this reason the position of Co G fio must be changed. We can 
think to correct its position and get a new point (cq + 5c) with a lower defect 
of invariance A = [1 — P]J(c + 5c). If the initial node is "not far" from the 
invariant manifold, a reasonable way to get the node correction 5c is to solve 
the linearized invariance equation where the vector field J is expanded to the 
first order and the projector P to the zeroth order: 

[1 - P(c)] [ J(c) + L{c)5c] = 0. (17) 

L is the matrix of first derivatives of J (Jacobian matrix). The Newton method 
with incomplete linearization consists of the equation (fT7j) supplied by the 
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extra condition [8]: 

P5c = 0. (18) 
The additional condition (j!8p and the atoms balances automatically can be 
taken into account choosing a basis in the subspace S = (kerP Piker D). 
Let h = dim(S), then the correction can be cast in the form 5c = Y$=i $ibi, 
so that the linearized invariance equation ([17]) becomes the linear algebraic 
system in terms of 8^. 

Ef=i 5 l ((1 - P)Lbi, b k ) = - ((1 - P) J, b k ), k = 1, fc. (19) 

Remark. Here the usual scalar product (, ) was used to get the components of 
the left-hand side of (|T7|) in the basis vectors Nevertheless, a different 

scalar product can be also used without a loss of generality. 

In the case of the thermodynamic projector, it proves convenient to choose 
the basis {bi} orthonormal with respect to the entropic scalar product (|HD and 
write the (IT3|) as: 

Eti 6t ((1 - P)Lb t , b k ) = - ((1 - P) J, b fc ) , fc = 1, fe. (20) 

The projector ffTB]) is "almost" (,) —orthogonal ((imP, ker P) = 0) close to 
the SIM. Because of that special feature, equation fliZO"]) can be approximated 
and simplified as follows: 

Eti <5< (Xbi, = - (J, b fc ) , k = 1, /i. (21) 

Note that, in general, an approximation carried out by eq. (121]) leaves a residual 
defect (1111) in the grid nodes which cannot be completely annihilated. There- 
fore, when a higher accuracy in the SIM description is required, equation (1191) 
is recommended. 



3.2.5 Iterative procedures: the relaxation method 



An alternative approach to solve eq. (fTTj) is the relaxation method. According 
to that method the correction is written as c = Co+r(c) A(c), and the quantity 
t(c) is obtained from the condition: 

(A,[l-P][J + r(c)LA]) = 0, 

and solving with respect to r: 

Equation (|22|) shows that the relaxation method is explicit, but as it adjusts 
the node position acting only along the direction of the defect A, typically we 
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expect it to be less efficient in comparison with the Newton method. On the 
other hand, this method is particularly easy to implement. 



4 The initial approximation. The Quasi Equilibrium Manifold 

Any iterative procedure needs to be supplied by a first approximation. Since 
it plays an important role for both the convergence and efficiency, that ap- 
proximation must be chosen very carefully. It was shown that a reasonable 
way, for initializing the MIG, is to construct the Quasi Equilibrium Manifold 
(QEM) [H2]. 

4-1 QEM definition 

Solution trajectories in the phase-space must obey the set of ODE equations 
([5]). Moreover, all trajectories also satisfy a subset of linear equations §6§ 
which represent the atom conservation. Among all the concentration points, 
that fulfill the latter constraints, we can choose those points which minimize 
the Lyapunov function G of the system we are dealing with. Such points 
lie on a manifold that is called Quasi Equilibrium Manifold (QEM). Let us 
suppose that some steps of a complex reaction are faster than some others. 
Since the Lyapunov function G must decrease during the fast dynamics then, 
when the fast motion is exhausted, the G value is expected to be the minimum 
on that fast hyperplane. In such a situation, a QEM attempts to achieve a 
motion decomposition into fast - toward the QEM - and slow - along the 
QEM. If the invariant manifold exists, the QEM can be taken as a reasonable 
approximation of it. In order to be more specific, let a chemical system have 
n reactive species. The degrees of freedom of that system are (n — I) because 
of the atom balances If q < (n — I) is the dimension of the QEM, then 
the macroscopic variables for its description are so that: (m b c) = 

£1, (m q , c) = £ g . Here, the n-dimensional vectors rrij are related to the 
hypothetic fast directions. From a mathematical standpoint, the solution of a 
variational problem: 



represents the QEM corresponding to the vector set {mj}. We want to stress 
the geometry behind (1231) because it will be extensively exploited in the fol- 
lowing. The geometric interpretation of a QEM is illustrated in Fig. [T] for a 




(23) 
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Fig. 1. Quasi Equilibrium Manifold: the geometrical interpretation. Two different 
QE- manifolds (bold lines in (a) and (b)) corresponding to two different linear con- 
straints subsets in the problem ([23]) . 



2-dimensional phase-space (c^., c^.)- Let us consider the points where G level 
curves (convex curves in Figures [TJ (a)-(b)) are cut by the QE-manifolds (bold 
curves): in those points the inclination of the tangent to the G-level curves 
is constant. Different QEM can be obtained by choosing different vector sets 
{rrii}. A special choice is done when m 1; ...,m q are the q left eigenvectors of 
the Jacobi matrix L(c eq ) corresponding to the q smallest absolute eigenval- 
ues. In that particular case, solution of ( |23|) has its own name: Spectral Quasi 
Equilibrium Manifold (SQEM) [T1I2] . 



4-2 Quasi Equilibrium Manifold in practice 



The minimization problem (1231) can be, in principle, solved by the method of 
Lagrange multipliers. However, it is also well known that, when the number of 
constraints and variables increases, then that method becomes prohibitively 
complicated to implement. Since the number of species and elementary reac- 
tion steps is usually quite high, the Lagrange multipliers method is not suitable 
for most of the practical cases. For this reason, in the sequel a new procedure 
to overcome that issue is presented. An algorithm (Quasi Equilibrium Grid 
Algorithm-QEGA) , which can be easily implemented in order to get a dis- 
crete analog of a QEM in any dimension, is developed. This is achieved by 
investigating further the QEM geometrical construction. 
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Fig. 2. Quasi Equilibrium Grid: the basic idea. 
5 ID Quasi Equilibrium Grid (QEG) construction 

Let us consider a one-dimensional quasi-equilibrium manifold. Let us assume 
that the node Cq belongs to that manifold. One may now imagine to look for 
a new node C\ which still lies on the quasi equilibrium manifold. In general, 
the node c± can be obtained from c by adding a shift 5c : c\ = c + Scq. 
That idea is applicable whenever a QEM-node c n is known and a new one 
c n+ i must be found (see Fig. [2]): 

c n +i = c n + 5c n . (24) 

First of all, any node c has to fulfill the atom balances (Q. Let {p^} be a basis 
in the null space of matrix D. A convenient way to take automatically into 
account the conditions ([6]) is to express any shift Sc n as a linear combination 
of vectors pf 

5c n = Y. Z i= l^P^ ( 25 ) 

where z = n — I is the dimension of the basis {pj}. By referring to Fig. 
[2J let us now discuss further the tangent space T to the G level surface in 
any quasi-equilibrium point c n+1 . The space T geometrically represents the 
linear constraint of the problem (I23I) . Therefore, any point c of T satisfies that 
constraint, but only c n+1 minimizes G function. The line f passing from c„ +1 
and c has the parametric form c = <y2t + c n+1 , where t is a vector of T spanning 
f while if is a parameter. In general, the linear constraints of the problem fl23l) 
can be also written as: 

(m, c) = <p(m, t) + (m, c„, +1 ) (m, t) = 0, Vt 

( 26 ) 

(dj, c) = <£>(di, t) + (dj, c n+1 ) ^ (dj, t) = 0, Vt 

where m and dj are the reduced variable vector (q = 1) and the generic row of 
matrix D, respectively. The vector t, that respects (12"oT) . can always be written 
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as a linear combination of some vectors tj, where {tj} denotes a basis in the 
null space of that matrix E, whose first row is given by m and the rest by the 
rows of D: 



E 



m 
D 



(27) 



Note that the dimension of {tj} is z — 1. By looking at Fig. [2J the quasi- 
equilibrium requirement simply becomes the orthogonality condition: 



(VG(c n+1 ),t) = 0, Vt GT (28) 

which also means: 



:vG( Cn+1 ),t,-; 



o, Vj 



2-1. 



(29) 



The quasi-equilibrium grid algorithm is based on the equation system (1291) 
and two more assumptions. First of all, we suppose that the known node c n 
is close to the QEM, although it does not necessarily belong to the QEM. 
Secondly, let the vector Sc n be small enough, so that the gradient VG(c n+ i) 
can be approximated to the first order: 



VG(ch-i) = VG(c n ) + H(c n )5c n 



(30) 



where H{c n ) = q c .q c . denotes again the matrix of second derivatives of the 
function G evaluated at the point c n . By substituting equations (1317]) and (1231) 
in (129]) . we obtain: 



EU (tj, H{c n ) Pi )in = -(tj, VG(c n )), Vj = 1, z - 1 • (31) 

By using the entropic scalar product <^j, equations (1311) can be cast into the 
form: 

E-=i (%, PillH = ~(tj, VG(c„,)), Vj = 1, z - 1 • (32) 

Both the matrix H and the gradient VG are calculated at the known node 
c n . Note that the right-hand side of (1321) vanishes if the node c n belongs to 
the QEM. The node collection, subsequently evaluated through (132]) . will be 
called a Quasi Equilibrium Grid (QEG). 



5.1 Closure through the spacing condition 



Note, however, that the system (1321) is not closed (z unknowns \x\, but z — 1 
equations) because it lacks a further information about the grid spacing. A 
reasonable closure for that system can be achieved by fixing the grid spacing 
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(e.g. in the Euclidean sense): 



5c n = e 



&,VGK)), Vj 



l,...,z 



(33) 



where e is a given number and 5c n represents the Euclidean norm of the 

vector 5c n . The smaller e is chosen, the more accurate the expression (1301) 
gets. As it will be shown later on, for small e the Quasi Equilibrium Grid 
lies very close to the correspondent Quasi Equilibrium Manifold. The extra 
condition makes (I33p a non-linear algebraic system. A way to solve it will be 
now discussed. The idea is to find the general solution of the linear system 
fl32l) . and then to choose the one which also fulfills the non linear condition in 
(1331) . Let the basis {p { } be orthonormal (in the Euclidean sense). That is not 
crucial, but it proves to be convenient in the following analysis; indeed the 
non-linear system (1331) now is cast as follows: 



(t i ,VC(c n )),Vj 



l,...,z 



(34) 



The general solution of (13*21) can always be written as: 



fix 




V\ 




Pi 




= w 




+ 












Pz 



"1, 



and p = [pi, 



(35) 



arc 



where w is a free parameter, while v 
the solution of the homogeneous problem and a special solution of (1521) . re- 
spectively. Without any restriction, we assume {y,v) = 1. Once f and p are 
known, the non linear equation of (|34l) can be written, in terms of w, as: 



w 



+ 2(is,p)w + (p,p) -e 1 



0. 



If the solvability condition is satisfied, 

(*/,p) 2 -(p,p)+ £ 2 >0, 



(36) 



(37) 



then the two real valued solutions of ( 1361) {w 1 , w 11 ), upon substitution into 



( 1351) . give two possible sets [//]_, //^]. Therefore, by using the ( 124]) and ( "251) . 
two new nodes c^+^c^ (both close to the quasi equilibrium curve) can be 
evaluated from the previous one c n (see Fig. [3]). A criterion, able to choose 
between those two solutions, depends on the phase-space zone where the grid 
needs to be constructed. This idea will be clarified in the following example. 
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Fig. 3. Two solutions for the ID QEG algorithm. 



Usually, the equilibrium point is supposed to be a good starting node for the 
QEG procedure: c = c eq . 

Remark. The QEG-equations (|32|) can be generalized as follows: 

Ef =1 (tj, Pi )m = -r)(tj, VG(c„)), Vj = 1, z - 1, (38) 

where r] is a parameter < rj < 1. When t] = 1, ( 1321) is recovered. On the other 
hand, if the QEG-nodes are close to the QEM, then the non-homogeneous 
terms can be neglected (they vanish on the QEM). Therefore, a reasonable 
approximation of the system (J32l) is given when rj — 0. In the latter case, 
the solvability condition ( |371) is fulfilled. If t] — 1 and ( |37j) does not hold, 
that parameter can be chosen in such a way that the solvability condition is 
satisfied. In the example below, solvability condition fl37|) is always satisfied 
and we use f l33|) . 



6 ID SQEG algorithm at work 



In this section, an example will be considered in order to illustrate how 
the algorithm, described in the previous section, works for finding a one- 
dimensional SQE-grid. That grid will be compared with the relative spec- 
tral quasi-equilibrium manifold, too. Let us consider the following four-step 
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Fig. 4. Reaction (i39j) : some trajectories projected into the phase-subspace (cc,cb)- 



three-component reaction (kindly suggested by A.N. Gorban): 



1. A++B, kf = 1, 

2. B «-> C, jfcj = 1, 
= 1, 

4.A + 5 «-> 2C, 4 



3.C «-> A, ^ 



(39) 



50. 



The atom balance takes the form: 



1 1 1 



ca 
cb 
cc 



(40) 



and the equilibrium point is chosen as: = 0.1, = 0.5, c§ = 0.4. As Fig. 

H] shows, the system is effectively two-dimensional, so the quasi-equilibrium 
manifold is expected to provide an one-dimensional reduced description (q = 

I) . Indeed, any solution trajectory, after a rapid initial dynamics, is attracted 
to a ID curve and along it reaches the equilibrium point. If the system is closed 
and the reaction fl39|) takes place under constant volume and temperature, we 
can assume that the system is supported by the Lyapunov function G ((7j): 



G = c A 




+ c B 




- 1 



+ c c 



In 



cc_ 



- 1 



(41) 
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Once a 3-dimensional vector m has been chosen, the QEM equation can be 
found by solving the variational problem (1231) : 



G — > min 
(m, c) = f 
Dc = 1. 



(42) 



In the following, the Spectral Quasi Equilibrium Manifold (SQEM) [2lfT] will 
be constructed. The Jacobian matrix L in the equilibrium point and its slowest 
left eigenvector xf are: 



L(c eq ) 



-30 -4.8 13.5 
-24 -6.2 13.75 
54 11 -27.25 



x. 



0.8807, -0.3905, 0.2681 



(43) 



Solution of the problem (|42l) . with the choice m = x s t , delivers the ID SQEM 
for the case shown in Fig HI To this end, let us rewrite the (H2l) in a more 
explicit form: 



c 0A = 0.3072 + 0.7867£ - 0.51800(0 
c 0B = 0.6928 - 0.7867£ - 0.48200(0 
coc = 0(0 



(44) 



gg(M) 



0. 



where c = [c a, Cob, c c] is the solution of the problem (l4"2"j) . while denotes 
the relation between cc and the reduced variable £ on the SQEM. By using 
the G function (jUJ), the problem (jUJ) is equivalent to the implicit equation 



'0.3072 + 0.7867^- 0.51800" 
OA , 



-0.5180 



'0.6928 - 0.7867£ - 0.48200 s 
05 . 



-0.4820 



0.4 



(45) 

The solution of fj4"5]) . by means of relations fT44l) . gives the SQEM shown in 
Fig. [5]^a). One may now apply the QEG-algorithm described above, in order 
to make a comparison with the analytic solution just found. An orthonormal 



basis {p^ in the null space of the matrix D 
and can be chosen as follows: 



1 1 1 



has dimension z 



Pl = [-0.5774,0.7887,-0.2113], 
p 2 = [-0.5774,-0.2113,0.7887]. 



(46) 
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0.3 0.5 0.7 0.9 1 



Fig. 5. (a) The bold curve is the SQE-manifold which was analytically evaluated 
by solving the (|45p . In that case the SQEM represents a very good approximation 
of the invariant manifold, (b) The SQE-manifold is compared with the SQE-grid 
where e 2 = 10~ 3 . 

Since the matrix E has the form: 



E 



0.8807 -0.3905 0.2680 
1 1 1 



(47) 



a vector t spanning ker(E) is: 



t = [-0.4229, -0.3934, 0.8163]. 



The system (1341) . in this example, simply reads: 



(t, pjm + it, p 2 )/i 2 = -(t,VG) 



(48) 



By solving (j4"8"|) in a QEG-node c n , the shift vector Sc n = H\P\ + \i2P2 allows 
to evaluate the new QEG-node c n+i = c n + 5c n . The QEG procedure, starting 
from the equilibrium point c eq = c , was performed twice, keeping uniformly 
parameter e 2 = 10" 3 . The first time, by choosing the solution in such a way 
that CB n+1 < cb„, the left branch of the SQE-grid was obtained; then, by 
imposing ce n+1 > ce n , also the right branch was calculated. The algorithm 
was terminated as soon as at least one component of the new node c n+ \ 
becomes negative. The result, shown in Fig. G^b), proves that the SQE-grid is 
in excellent agreement with the analytical curve (SQEM). 
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Fig. 6. SQEG Left branch of the case in Fig. [5] (b). Different approximations com- 
pared with the analytical solution (SQEM). Each grid is calculated by using a 
different parameter e. 



6.1 Grid spacing choice 



There is no need to stress the importance of the grid spacing parameter e for 
the QEG accuracy. In the case of Section [6J the SQEG was computed several 
times with different values of e. The QEG algorithm is based on the linear 
approximation (1301) . Therefore, the smaller is 5c n = e the more accurate is 
the QEM description by means of the QEG. Nevertheless, the smaller is e the 
larger is the number of times that the system fl33|) must be solved to have a 
grid of a fixed size. For this reason, we need to keep e as large as possible. 
We estimated (at least the order of magnitude) the upper limit of spacing 
(e u ) which gives a QEG "not far" from the relative QEM. From our numerical 
experiments, a reasonable value for that was e u = 10 _1 . As Fig. [6] shows, the 
QEG is not far from the QEM even for a quite coarse grid (e > e u ). 



7 Generalization to multi-dimensional grids 



The QEG algorithm, which has been developed for constructing 1-dimensional 
grids, can be modified in order to get multi-dimensional grids, whenever 
needed. From all reasonable extension strategies, two of them here will be 
analyzed: a straightforward extension and a flag extension (flag extension, for 
invariant grids, was introduced in Ref. [3]). In the first case, the algorithm 
of paragraph and the equation system (|34[) are tuned for a g-dimensional 
grid calculation. Here, the implicit assumption is that the grid dimension q is 
fixed and uniform everywhere in the phase space (like for the QEM construc- 
tion). However, a second flexible approach, suitable for SQEG construction, 
was developed, too. In that case, the grid dimension can be varied at will. 
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7.1 The straightforward extension 



According to the straightforward extension, if a node c n close to the q- dimensional 
QEM is known, then a new node c n+1 can be added to the QE-grid by shifting 

Cn+i = c n + 5c n , 5c n = ^2 i=l fiiPi, (49) 

where {p^} is still a basis in the null space of matrix D. The linear constraints 
of the problem (1231) define the tangent space T to the G level surfaces in the 
new node c n+1 . Let c be a generic point of T, the line f passing from c n+1 and 
c has the parametric form: c = tpt + c n+1 , where t is the vector of T which 
spans f and 99 is the parameter. The generalized form of the relations (126]) is: 



mi, c) =(p [mi, tj + (mi, c l n+1 j [m 1 , tj = 0, Vt G T 
m 9 , c) = y? (m„ tj + (m q , =4> (m q , tj = 0, Vt G T 

(di, c) = <f (d h tj + (dj, ^ (di, t) = 0, Vt e T, 



(50) 



which means that vector t belongs to the null space of the matrix E {keiE)\ 



E 



mi 



D 



(51) 



Now, the dimension of basis {tj} in ker(.E) is (z — q). Since the quasi equi- 
librium condition requires that, among all the points c of T, c n+i has the 
minimal value of G, the following orthogonality conditions hold: 



;VGM,f 3 ) = 0, Vj = l,...,z-q. 



(52) 



For small vector 5c n , the approximation (1311 can be used, so that the (1521 
become: 



£* =1 Pi) IH = - VG(c„)) , Vj = 1, « - g. 



(53) 



As the system (15B"|) shows, the larger is the QEM dimension (5) the smaller 
is the set of "mere" quasi-equilibrium equations available, while the number 
of unknowns remains constant (z). The closure of the rectangular system (1531 
requires q more equations and has only to do with the geometric structure 
which we want to provide the grid with (e.g. grid spacing, shift vector ori- 
entation in the phase-space, etc). In general, the geometric structure of the 
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2D Quasi Equilibrium Manifold 



Fig. 7. 2D Quasi Equilibrium Manifold. Location of solutions of the system (|54p in 
the phase space. 

grid under construction can be chosen at will: therefore there is no unique 
geometric closure for that system. However, one possible condition could be 
imposed, like in fl33l) . by fixing the Euclidean norm of shift vector: 5c n = e. 
Nevertheless, (q — 1) geometric constraints are still missing. In order to illus- 
trate how the geometric closure issue can be overcome, the case q = 2 will 
be considered in the following. For that special possible closure, which 

can be easily generalized, will be presented. If a two-dimensional QEG has to 
be constructed, then only one extra equation is needed to close the system: 



EU Pi) fM = - (tj, VG(c„)) , Vj = 1, ...,z - 2 

5c r 



(54) 



Fig. [7] shows that all the possible solutions of (154"1) are located, as a "crown" , 
near the QEM. A way to choose only two of them can be achieved by intro- 
ducing a new fixed vector rh and imposing a given angle between rh and 
5c„: 



J2 i= i (™> pi) ^ 



5c r 



\rh\\ cos$. 



(55) 



The choice $ = n/2 proves to be particularly convenient, as ( 1551) becomes: 

5D* = i Pi) Vi = °- ( 56 ) 

( 1561) allows to write a closed system: 

Ef=l (tj, Pi) IH = - (tj, VG(Cn)) , Vj = 1, ...,z - 2 



(57) 



<5c r 



where the extra information, through e and rh, concerns the grid spacing 
and the phase-space zone of interest where the grid must be constructed. In 
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general, the geometric closure of (155]) can be achieved when (q—1) independent 
vectors {m-j} and the parameter e are fixed. Here, we present an approach 
which allows to get a rectangular structured grid. The general form of (157)) is: 

EiU Pi) ^ = - VG(c„)) , Vj = 1, * - q 
< Eti (Aj, Pi) IH = 0, Vj = 1, q - 1 (58) 
5c n = e. 

The g-dimensional grid construction is split in q subsequent steps. Starting 
from the equilibrium c eq , system (1581) is solved by choosing (q—1) rrij vectors 
among the q available and imposing: rhj = rrij Vj = l,...,q — 1. In this 
way, a first set of QEG nodes is attained as soon as e is known. Now, starting 
from each of those points, system (1551) . by using a different combination of rrij 
vectors, gives some more nodes. The procedure ends (g-th step) when all the 
possible different combinations of (q — 1) vectors {rrij} are over. In Section [HJ 
that idea will be explained by means of an illustrative example. 

7.2 The flag extension 

A multi-dimensional QE-grid construction becomes non-trivial especially when 
q becomes large. As reported in section 17.1) the straightforward extension re- 
quires to introduce some additional vectors rhj. The flag extension can be 
applied when a SQEG is searched. That procedure is strongly based on the 
algorithm presented in paragraph [5J and it naturally leads to a rectangular 
structured grid. The idea, which is behind, is simple and makes this method 
very flexible and really suitable for constructing high-dimensional rectangu- 
lar grids. Let us suppose that q is the grid dimension and the q SQE-vectors 
{mi, ...,m s } are fixed. Here, the assumption is that m 1 is the slowest eigen- 
vector (corresponding to the smallest eigenvalue by absolute value), m<i the 
second slowest and so forth. The grid construction is achieved in s subse- 
quent steps. In each step one more dimension is added to the grid. At the 
beginning, by using m = mi, the algorithm in section [5J provides the ID 
quasi-equilibrium grid. Now, starting from any node c* of that grid, a new ID 
QEG is constructed where m = rri2. In this case, the second QEG represents 
a trajectory on the 2D-manifold attracted to the slowest ID-manifold in the 
node c*, once the fast dynamics is exhausted (see Fig. [8]). G function depends 
on the equilibrium point c eq : G = G(c,c eq ). Since c* can be considered as 
a "local equilibrium" for the fast motion, the second ID grid is obtained by 
minimizing G = G(c,c*). Once the previous step is completed, the grid can 
be extended in the third dimension by adding, in each node c of the new 2D 
grid, a ID QEG where m = m 3 and G = G(c, c ). In this way, the procedure 
is performed up to a g-dimensional grid. By extending partially a previous 
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Fig. 8. A 2D flag. Once the ID quasi equilibrium grid is found, from each node c*, 
new ID quasi equilibrium grids are added. The second slowest ID grid represents 
that trajectory collected by the first ID quasi equilibrium grid in the node c*. 

grid, it becomes really easy to have some grids whose dimension is different 
in different phase space zones. It is worth to stress that the straightforward 
and the flag extension deliver two different objects: the first one just gives 
the quasi-equilibrium grid "brute force", while the second one is its conve- 
nient "approximation" which has some useful features as it will be illustrated 
in the following. First of all, the flag grid does not demand any extra vec- 
tor for the geometric closure and the grid dimension can be easily varied in 
different phase-space zones. Secondly, if a grid refinement procedure (MIG) 
is used in order to get an invariant grid out of the quasi-equilibrium one [2], 
then the flag extension reveals to be a very useful tool. Indeed, let us assume 
that a multi-dimensional invariant grid is required in order to reduce a given 
model. A possible strategy might be given by a "hybrid procedure" where the 
QEG algorithm and the MIG method are alternatively used according to the 
sequence: 

• ID quasi-equilibrium grid construction (slowest grid); 

• MIG refinements until the ID invariant grid is obtained; 

• flag extension from ID invariant grid to 2D quasi-equilibrium grid; 

• MIG refinements until the 2D invariant grid is obtained; 

• flag extension from 2D invariant grid to 3D quasi-equilibrium grid; 

• MIG refinements... 



7.3 Beyond SQEG: GQEG and SEGQEG 

The latter suggestion sheds light on one more option which, if implemented 
during the flag extension, allows to go beyond the SQEG approximation of the 
invariant manifold. Let us assume that the hybrid procedure of Section 17^21 is 
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utilized and a k- dimensional invariant grid (let c* be its generic node) has to be 
extended to a (A; + 1) -dimensional grid. That grid will approximate the (k + 1)- 
dimensional invariant grid, better than the SQEG does, if in each invariant 
node c* the vector m is chosen as the (k + l)-th slowest left eigenvector (by 
absolute value) of Jacobi matrix L(c*). According to [Tj, here a considerable 
simplification can be achieved by replacing the full Jacobian L(c*) with: 



too. Matrix L sym is symmetric with respect to the entropic scalar product 
dSJ). For that reason the spectral decomposition will be much more viable (see 
also Ref. [2]). Those two new approximations will be named: Guided Quasi 
Equilibrium Grid (GQEG) when the full Jacobian L(c*) is used, while Sym- 
metric Entropic Guided Quasi Equilibrium Grid (SEGQEG) if L sym replaces 
the full matrix. In order to give an idea about the effort needed, for example 
in a SEGQEG construction, let us consider a 2-dimensional grid. In that case, 
the spectral decomposition of a symmetric operator is performed only over 
the nodes of an one-dimensional grid. Moreover, also a possible criterion, for 
getting a multi-dimensional grid, naturally applies: if at the node c* of the k- 
dimensional invariant grid, the ratio |Ajt+i|/|Afe| (between eigenvalues of L or 
X sym , respectively) is not larger than a fixed threshold, the (/c + l)-dimensional 
grid will not be extended at that point. In this way, the grid dimension q is 
generally not uniform in the phase space. Several techniques suggested above 
are only some reasonable ones. The flexibility of the method proposed allows 
to set up different procedures, still based on the Quasi Equilibrium Grid ap- 
proach: the QEG system fl53|) supplied by a geometrical closure. 



8 2D Grid Example: hydrogen oxidation reaction 



Let us consider a model for hydrogen oxidation reaction where six species H 2 
(hydrogen), 2 (oxygen), H 2 (water), H, O, OH (radicals) are involved in 
six steps in a closed system under constant volume and temperature (see Ref. 




(59) 



where L T is the ordinary transposition, and H is evaluated at the point c*, 
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IU, P- 291): 



l.F 2 «-> 2F, kt = 2, 



2.0 2 «-> 20, fcj = 1, 



c 3~ — 1) 



3. F 2 «-> F + OF, ^ 

4. F 2 + O <-> F + OF, jtf = 10 3 , 
5.0 2 + F <-> O + OF, = 10 3 , 
6.F 2 + O <-> F 2 0, fct = 10 2 . 



(60) 



The conservation laws are: 



2ch 2 + 2ch 2 o + c H + c H = b H = 2 
2c 02 + c H2 o + c + c h = b Q = 1. 



(61) 



When the equilibrium point is fixed, for example 



eq _ n 2? eq 
C H 2 — U.Z/, C 02 



0.135, c% = 0.7, 4 9 = 0.05, c% = 0.02, c e J H 



0.01, 



(62) 

then the rest of the rate constants k i are calculated using the detailed balance 
principle (jlj). The system under consideration is fictitious in the sense that 
the subset of equations corresponds to the simplified picture of this chemical 
process and the rate constants reflect only orders of magnitude for relevant 
real- word systems. We can assume that the Lyapunov function G has the form: 



g = E 6 1 



In 



- 1 



(63) 



Here, we are interested in the 2D SQEG construction. Two left eigenvectors 
of Jacobian matrix L(c eq ) are: 



xf = [-0.577,-0.568,0.225,0.0482,0.0666,-0.536] 
xf = [0.00682, -0.00595, 0.0221, -0.7, -0.713, 0.423] 



(64) 



where xf 1 and xf 2 are the slowest and the second slowest one, respectively. 
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Fig. 9. The 2D SQEG constructed by using the straightforward extension with 
e 2 = 0.5 • 10~ 3 : projection into the phase-subspace (ch,cq,coh)- 



8.1 The 2D straightforward extension 



In order to get a 2D SQEG for that example, the straightforward extension 
was used as first strategy. Matrices D and E take now the form: 



-0.577 -0.568 0.225 0.0482 -0.0666 -0.536 

0.00682 -0.00595 0.0221 -0.7 -0.713 0.423 

2 2 1 1 

2 10 11 

(65)" 

As suggested in the end of section [7TT| the procedure has been started from the 
equilibrium point and it was split in two subsequent steps. At the beginning, 
the system (1571) was solved by imposing e 2 = 0.5 • 10~ 3 and rh = xf 2 : in this 
way, the grid nodes, denoted by circles, in Fig. [9] were obtained. In the second 
step, (1571) was solved by starting from any circle: this time, the geometric 
constraints were e 2 = 0.5 • 10~ 3 and rh = xf 1 . During that step, in each circle, 
the horizontal dots of Fig. [9] were found, too. 



D 



2 2 10 1 
2 10 11 



E 
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Fig. 10. ID Spectral Quasi Equilibrium Grid with e 2 = 3 • 10 3 : comparison with 
the ID invariant grid obtained by MIG refinements. 

8.2 The 2D flag extension 



After that, also the flag extension procedure was applied. Now, the ID spectral 
quasi equilibrium grid is needed. Matrices D and E are in this case: 



D 



2 2 10 1 
2 10 11 



E 



-0.577 -0.568 0.225 0.0482 0.0666 -0.536 
2 2 1 1 

2 10 1 1 



(66) 

Starting from the equilibrium point c eq , the system (134"]) was solved by fixing 
e 2 = 3 • 10~ 3 (see Fig. [TU]). The flag extension was used to get a 2D grid out 
of the ID one. Now, the new matrix E reads: 



E 



0.00682 -0.00595 0.0221 -0.7 -0.713 0.423 
2 2 1 1 

2 10 11 



while the Lyapunov function G has the form: 



G 



E 6 



hi 



Cj 



- 1 



(67) 



where c* = [c*,...,Cg] is any ID grid node which is extended in the second 
dimension (see Fig. [H]). Figures [TITa)-(b) show two different 2D SQE-grids: 
the first one is obtained by extending the ID SQE-grid, while in the second 
case the ID invariant grid is used. In other words, the latter result was attained 
by the "hybrid procedure" QEGA + MIG suggested in the end of section I7T21 
For both cases, in the second dimension, the grid spacing was e 2 = 1.5 • 10~ 3 . 
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Fig. 11. The flag extension, (a) 2D SQEG (dots) extended from the ID SQEG 
(circles), (b) 2D SQEG (dots) extended from the ID invariant grid (circles). The 
grid spacing, in the second dimension, was e 2 = 1.5 • 10 -3 . 

8.3 The 2D GQEG and SEGQEG 

Finally, the GQEG and SEGQEG approximations are computed for the hy- 
drogen oxidation reaction ( 160]) (Lease). Here, the grid spacing is uniformly 
kept e 2 = 0.45 • 10~ 3 . Each grid has 10 x 15 nodes and it is compared with 
both the SQEG (straightforward extension) of similar size (check Table [1]) and 
the invariant one. The invariant grid was obtained by refining the approxima- 
tions through the MIG procedure. All those grids lie quite close to each other. 
However, a "more pathological" case (2. case) is also analyzed (here the SQEG, 
far from the equilibrium, presents a remarkable deviation from the invariant 
grid): now the rate constant set is taken as kf — 20, k£ = 1, k£ = 1, k£ — 10 3 , 
k§ = 10 3 , k$ = 10 2 , while the equilibrium point coordinates still are given 
by (J62"j) . For that case, the SQEG, GQEG and SEGQEG were constructed by 
choosing the grid spacing and size as for the previous case. Note that all the 
grids were partially extended only below the equilibrium point in the phase- 
space zone where they present the largest deviation from the invariant grid. 
This time those three approximations have a low invariance defect only near 
the equilibrium. In order to estimate how far each grid is from the invariant 
one, the following procedure is implemented. A 10 x 15 matrix, collecting in 
any grid node an invariance defect measure, is constructed. As suggested by 
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002 0.03 0.04 



0.06 0.07 



0.08 0.09 0.1 



Fig. 12. Lease: k\ = 2, k% = 1, k£ = I, k£ = 10 3 , k$ = 10 3 , fc^ = 10 2 , 
e 2 = 0.45 • 10 3 . Two grids formed by 10 x 15 nodes. A 2D GQEG (dots) and a 
2D invariant grid (circles) are reported. The grids are partially extended below the 
equilibrium point (square). 




O 0.02 0.04 0.06 0.08 0.1 .12 0. 



14 0.16 0.18 



Fig. 13. 2.case: kf = 20, jfej 



1, kn 



1, k/\ 



io 3 , kt = io 3 , k: 



io 2 



e 2 = 0.45 • 10 3 . Two grids formed by 10 x 15 nodes. A 2D GQEG (dots) and a 
2D invariant grid (circles) are reported. The grids are partially extended below the 
equilibrium point (square). 



[2], that local measure may be y (A, A) / (J, J), where A and J are the in- 
variance defect ([PIT) and the vector field of (jSJ), respectively. A is evaluated by 
using the thermodynamic projector (jl~6]) . By averaging over all the invariance 
defect measures, the mean invariance defect is provided: results for both cases 
are condensed in Table [H Note that the adopted invariance defect measure is 
dimensionless as it compares the invariance defect with the vector field. Cal- 
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1 case 


2 case 


SQEG 


0.318 


0.645 


GQEG 


0.238 


0.460 


SEGQEG 


0.303 


0.491 



Tablp 1 

Mpan invariancp dpfpct (dimpnsionlpss): thrpp approximations of the invariant grid 
under comparison for the hydrogen oxidation reaction. In Lease, the parameter set 
is: k+ = 2,k+ = 1, fc+ = 1, fc+ = 10 3 , k+ = 10 3 , k+ = 10 2 , e 2 = 0.45 -10 3 . In 2. case, 
the parameter set is: = 20, k% = 1, k£ = 1, = 10 3 , £^ = 10 3 , /cj = 10 2 , 
e 2 = 0.45 • 10 3 . 

culations prove that the GQEG is better than the SQEG (straightforwardly 
extended); nevertheless the SEGQEG construction, since it requires a much 
lower computational effort and still has an error similar to the GQEG, is 
recommended when the SQEG is considered not satisfactory (e.g. big mean 
defect). 



9 Conclusions 

In this paper, the problem of Quasi Equilibrium Manifold approximation by 
means of a grid description is addressed. To this end, the notion of Quasi Equi- 
librium Grid (QEG) is introduced and a proper algorithm to construct it, in 
any dimension, is suggested (QEGA). It has been shown, through illustrative 
examples, that the QEGA gives a very good QEM approximation without fac- 
ing the analytical difficulties of Lagrange multipliers method implementation 
in large dimension. Since the QEGA is a completely numerical procedure, it 
reveals particularly suitable for providing the MIG procedure with the first 
SIM approximation. As it has been illustrated, some proper hybrid procedures 
QEGA + MIG, where both methods are alternatively used, allow to obtain 
accurate SIM approximations. It was proved that two special QEGA + MIG 
procedures deliver enhanced approximations of SIM: the Guided Quasi Equi- 
librium Grid and the Symmetric Entropic Guided Quasi Equilibrium Grid. 
Here, we want to stress the two major advantages of the method proposed. 
First of all, it is a completely numerical algorithm which only deals with nodes 
sets. Moreover, it is a local construction: namely, the computation of a new 
node Cn+i, which has to be added to the grid, only depends on the previous 
neighbor c n . Those two points make the QEG construction suitable for numer- 
ical applications and parallel realizations. Finally, it is not excluded that the 
QEGA is applicable not only for model reduction, but in some very different 
fields, too. Indeed, it was mentioned that the QEM notion already is exploited 
for some applications in Lattice Boltzmann schemes simulations. More gen- 
erally, the QEGA is a numerical tool which can be used to find a grid-based 
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approximation for the locus of minima of a convex function under some lin- 
ear constraints. In this paper we focused only on the geometry of the model 
reduction, that is, construction of slow invariant manifolds approximations. 
The implementation of grid-based integrators for dynamic equations will be 
presented in a separate publication. 
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